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It is widely accepted that both ripples and dunes form in rivers by primary linear in- 
stability the wavelength of the former scaling on the grain size, that of the latter being 
controled by the water depth. We revisit here this problem, using the computation of 
the turbulent flow over a wavy bottom performed in Part 1. A multi-scale description 
of the problem is proposed, in which the details of the different mechanisms controlling 
sediment transport are encoded into three quantities: the saturated flux, the saturation 
length and the threshold shear stress. Theses quantities are modelled in the case of 
erosion and momentum limited bed loads. This framework allows to give a clear picture 
of the instability in terms of dynamical mechanisms. The relation between the wave- 
length at which ripples form and the flux saturation length is quantitatively derived. 
This solve the discrepancy between measured initial wavelengths and predictions that do 
not take this lag between flow velocity and sediment transport into account. Inverting 
the problem, experimental data is used to determine the saturation length as a function 
of grain size and shear velocity. Finally, using the systematic expansion of the flow field 
with respect to the corrugation amplitude performed in part 1 , we discuss the non-linear 
selection of ripple aspect ratio. Investigating the effect of a free surface on the linear 
instability, we show that the excitation of standing waves at the surface has a stabilising 
effect, independently of the details of the flow and sediment transport models. Conse- 
quently, the shape of the dispersion relation obtained from the linear stability analysis 
of a flat sand bed is such that dunes can not result from a primary linear instability. We 
present the results of field experiments performed in the natural sandy Leyre river, which 
evidence the formation of ripples by a linear instability and the formation of dunes by 
a non-linear pattern coarsening limited by the free surface. Finally, we show that mega- 
dunes form when the sand bed presents heterogeneities such as a wide distribution of 
grain sizes. 



1. Introduction 

Since IRichards 19801 and followed by others QSumer fc Bakioglu 1984} IMcLean 1990| , 
ripples and dunes observed on the bed of sandy rivers have been interpreted as the 
two most unstable modes of the same linear instability. Although the classification of 
river bedform types is a difficult task ( |Ashley 1990| ), subaqueous ripples and dunes are 
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standardly defined by their typical size: ripples are the small scale bedforms whose 
wavelength A scales on the grain size and dunes are those whose wavelength is compa- 
rable to or larger than the flow depth H ( Kennedy 19 63 Richards 1980, Engcl und 1970[ 
|Freds0e 1974[ [ Engelund fc Freds0e 1982} lAllen 1985[) . A variant of this definition was 
used by |Guy et al. 1966 dunes are the bed features larger than ripples that distort sig- 
nificantly the free surface and are out of phase with the standing waves they generate; rip- 
ples are the triangular shaped bed-forms that have "lengths of less than about 2 feet and 
heights of less than about 0.2 foot". In some other articles (Hill et al. 1969, Yalin 1977), 
the classification of bedforms is rather based on the ratio of the grain diameter d to the 
viscous sub-layer depth, i.e. on the particle Reynolds number u^djv (u* is the shear 
velocity and v the kinematic viscosity of the fluid). These criterions are questionable, 
especially when there is no clear separation of length-scale between the spacing of ripple 
crests and the flow depth, and also in the case of bedforms superimposed on larger struc- 
tures (jVenditti et al. 200 5a). In this paper, willing to refer to the dynamical mechanisms 
involved in the formation of these bedforms, we shall call ripples those whose charac- 
teristics are independent of the flow depth; in contrast, dunes are directly related to 
the presence of a free surface. In particular, we shall not use the terms 'sand wavelets', 
as suggested by Coleman et al. (jColeman fc Melville 19941 IColeman fc Melvi lle 1996, 
|Coleman fc E ling 2000] IColeman et al. 2003|l to designate the initial stage of ripples 
when they emerge from a flat sand bed. 

Experimental measurements on subaqueous ripples are numerous and exhibit large 
data dispersion. One of the reasons for such a dispersion is that ripples exhibit pattern 
coarsening, i.e. present a progressive increase of their typical length-scale as time goes by 
(jMantz 1978l|Gyr fc Schmid 1989|lBaas 1994IIColeman fc Melville 1994IIColeman fc Melville 19961 
IBaas 19991IRobert fc Uhlman 200TllColeman et al. 2003UVenditti et al. 2005bl [Valance fc Langlois 2008} 
IRauen et al. 2 008). Many authors have only measured fully-developed wavelengths, 
whose relation with the initial wavelength A is still an open issue that will be discussed 
here. A reference data base of final wavelengths in inclined flumes has been completed 
by lYalin 19851 Although not directly comparable to the predictions of a linear stability 
analysis, these data show distinct scaling laws for hydraulically smooth and rough gran- 
ular beds. For a particle Reynolds number u*d/v smaller than few unities, the viscous 
sub-layer is larger than the grain size d and the fully-developed ripple wavelength turns 
out to scale on v/u*, with a large - yet unexplained - prefactor (~ 10 3 ). On the other 
hand, when the grain diameter d is larger than the viscous sub-layer, the wavelength 
was found to slowly increase with it*. Here we will focus on this hydraulically rough 
regime, in which the flow is turbulent at all scales down to the grain size. This cor- 
responds to the domain of validity of the calculations performed in part I. In several 
experimental articles (|Baas 1994| IColeman fc Melville 19941 IColeman fc Melville 19961 
|Coleman fc Eling 2000[ IColeman" et al. 20031 [Valance fc Langlois 20051 , the initial wave- 
length A has been carefully determined, showing scaling laws independent of the flow 
depth H . At large particle Reynolds number, A is found to be almost independent of the 
shear velocity u* and to increase with the grain diameter d (for instance lOoleman et al. 20031 
propose A oc d 75 ). These measurements will serve as a benchmark for the theory devel- 
oped in the present paper. 

Concerning dunes, most of the measurements have been performed in natural rivers, 
for which the Froude number T is low (see IBest 20051 for a recent review). In the 
Mississipi river ([Harbor 1998|) . in the Missouri river (|Ann ambhotl a"T972[l . in the Rhine 
dCarling et al. 2000[IWilbers fc Ten Brinke 2003[) and in the Rio Parana (|Parsons et al. 2005|) . 
the observed wavelengths range from 0.57? to 20-ff for Froude numbers around 0.2. The 
suggestion by |Kennedy 1963[ lAllen 19851 or IColombini 20041 that mature dunes should 




Figure 1. Schematic of the instability mechanism showing the streamlines around a bump. 
The fluid is flowing from left to right. A bump grows when the point at which the sand flux 
is maximum is shifted upwind the crest. The upwind shift of the maximum shear stress with 
respect to the crest scales on the size of the bump. The spatial lag between the shear and flux 
maxima is the saturation length L Ba ,%. 



present a well selected wavelength scaling on HjT 2 is thus far from reflecting the natural 
dispersion of field data. Flume experiments (see |Guy et al. 1966 IRobert fc Uhlman 20011 



and the data collected by |Kennedy 1963 1 have been performed in a much larger range 



of Froude numbers (from 0.1 to 1) and also show well-developped bedform wavelengths 
between H and 30_ff . The extensive data published by |Guy et al. 1966| are particularly 
impressive, and the interested reader should refer to the original article rather to the 
truncated data-set plotted by |Freds0e 1974| We will discuss these data in detail in the 
last section of the present article. 

It is widely accepted that ripples result from a linear instability whose key mechanism 
is the phase difference between sediment transport and bed topography ( |Kennedy 1963[ 
Reynolds 1965l|Kennedy 1969ilSmith"T970l|Hayashi 1970|lPa7ker 1975HEngelund fc Freds0e 1982| 



McLean 19901) . This spatial shift can be conceptually splittcd into two contributions. As 



seen in Part 1, there is a lag between the bed elevation profile and the basal shear 
stress, encoded in the ratio B/A, which results from the simultaneous effects of inertia 
and dissipation i.e. from hydrodynamics only (Fig. [T]). In addition, sediment trans- 
port needs some time/length to adapt to some imposed shearing. Kennedy 1963) and 



his followers ( |Engelund 1970| |Hayashi 19701 ISmith 19701 |Freds0e 19741 ) have introduced 
a phenomenological spatial lag L sat (this lag is called i5 in these papers) in a phase shift 
equation: 

q(x) = q s& t(x - L sat ), (1.1) 

where g sa t is the saturated (or equilibrium) flux, which is a function of the basal shear 
stress. With such relation between q and q sa t, the phase shift between the sediment 
flux and the shear stress over a relief of wavelength A is then simply — 27rL sat /A, which 
corresponds to a phase delay for A > 2L sat , but to a phase advance for L sat < A < 2L sat . 
As shown bv lParker 19751 such a constant phase lag is not physically founded and should 
be replaced by a first order relaxation equation of the form 

L sat d x q = <7sat - q- (1.2) 



This linear equation reflects the fact that the sand flux reaches its saturated value over 
a characteric length L sat ([Parker 19751 ISauermann et al. 20011 lAndreotti et al. 2002a[ 
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lAndreotti et al. 2002bl|Kroy et al. 2002||Wia"nce fc Langlois 2005| [Ch arru 2006|) . As for 
any other first order system, the phase delay between the flux and the shear stress 
increases here from for large wavelengths to 7r/2 for wavelengths much smaller than L sat . 
When the sum of the phase lags coming from hydrodynamics and transport results into a 
maximum flux upstream the bed crest, sediment deposition occurs on the bump, leading 
to an unstable situation, i.e. to the amplification of the bedform (Fig. [T]). Conversely, 
when the maximum flux is located downstream the crest, the bump is eroded and the 



disturbance decay. With the flux equation ( 1 . 1 1 , even a potential description of the flow 
( |Kennedy 1963[ |Reynolds 1965[ ) for which the basal shear stress is in phase with the 
topography (see part 1), can then find an unstable - but unphysical - region. However, 
sediment transport has a definite stabilising role, and the destabilising mechanism must 
result from hydrodynamics. 

The description of a turbulent flow over a sinusoidal bottom has been, in this context, 
progressively improved by |Engelund 197O}|Freds0e 1974[ [Richards 1980l and |Sumer fe B akioglu 1984. 
None of these articles take into account a saturation length. As discussed in Part 1, the 
basal shear stress is generically found to be in phase advance with respect to the bot- 
tom profile. This means that, without any stabilising mechanism, all wavelengths are 
found unstable. A key mechanism introduced by |Freds0e 1974| and kept in the follow- 
ing models is the stabilisation of short wavelengths by slope effects. In these models, 
the prefered wavelength results from hydrodynamics and is proportional to the hy- 
drodynamical roughness zq (or to vju* in the case of a viscous surface layer treated 
by |Sumer fc B akioglu 1984]). As shown bv ICharru 2006[ the predicted wavelengths are 
smaller than all experimental findings by several orders of magnitude. 

The concept of flux saturation was correctly introduced in a linear stability analysis by 
Parker 1975, the flow being modelled by depth averaged equations standardly used in hy- 
draulics. However, these equations cannot describe the layered structure of the turbulent 
flow above bumps (see [Jackson fc Hunt 19751 and Part 1), and in particular the so-called 
inner layer which is responsible for the phase shift between shear and topography. As a 
consequence, this author missed the explaination of the ripple instability. The idea that 
-L sat is the relevant lengthscale for the problem of dune formation was introduced in the 



aeolian context bvlSaucrma nn et al. 200ll lAndreotti et al.~2 002b Kroy et al. 2002 



ear stability analyis developed in this context (Elbelrhiti et al. 2005, IClaudin fc AndreottT 2006 



Lin- 



lAndreotti fc Claudin 200"7| suggest that the initial wavelength actually scales on L sat , 
and not on z . This line of thinking has been applied since then to subaqueous rip- 
ples (jValance 20051 |Valance fc Langlois 2005[ ICharru 2006[) . Although L sat and zq are 
both ultimately related to the grain size d, they correspond to very different mechanisms 
associated to sand transport and hydrodynamics respectively. 

Although Richards 1980 suggested that dunes could form by the same linear instability 
as ripples, no article has ever exhibited a proper and complete dispersion relation showing, 
in the same graph, the growth rate a for a range of wave-numbers fc that includes 
both ripples and dunes. In the figures provided bv Richards 1 980 a is rescaled by 
fc 2 , which artificially enhances the small fc (dunes) region. More recently, Colombini 
(Colombin i 2004[ [Colombini & Stocchino 2005) revisited this problem, introducing the 
thickness of the transport layer has a key parameter. This empirical refinement gives in 
practice an adjustable phase shift between relief and sand transport and allows to get a 
single well-defined peak in the dispersion relation, associated to dunes, at a wavelength 
on the order of the flow depth. However, this model does not predict the existence of a 
small-scale (ripple) instability anymore. 

In this paper, we wish to discuss afresh this ripple and dune formation problem, and in 
particular to get such a complete relation dispersion, where the peak associated to ripples 
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is present. To reach this goal, we must first discuss the sediment transport issue, before 
mixing it with the hydrodynamics presented in Part 1. In contrast to IRichards"! 980 
and IColombini 20041 we eventually show that dunes cannot form by a linear instability 
mechanism, and must rather result from the interaction and non-linear pattern coarsen- 
ing of ripples, limited by the free surface, as suggested by Raudkivi & Witt e 19901 and 
IRaudkivi 20061 

This paper is structured as follows. In the next section devoted to the sediment trans- 
port issue, we show how to abstract the details of the transport mechanisms into a general 
framework for the description of the different modes of transport and thus valid in the 
subaqueous and aeolian cases. A quantitative model of the transport threshold that in- 
cludes the slope effects is proposed. We then show that the sediment transport can be 
limited either by the erosion rate or by the momentum available in the flow. With these 
ingredients, we propose a self-consistent model that unifies these two mechanisms and 
discuss the scaling law of the saturation length L sat . In section 3, we revisit the forma- 
tion of current ripples (and aeolian dunes) by linear instability, using the quantitative 
hydrodynamic calculation performed in part 1 . It confirms that the most unstable wave- 
length is determined by the saturation length, and not by the hydrodynamical roughness. 
Inverting the problem, this allows to determine the saturation length from experimental 
measurements of the initial wavelength and to discuss its dependence on the grain size and 
the shear velocity. Finally, using the systematic expansion of the flow field with respect 
to the corrugation amplitude performed in part 1, we discuss the non- linear selection of 
ripple aspect ratio. In section 4, we investigate the effect of a free surface on the linear 
instability. We show that the excitation of standing waves at the surface has a stabilising 
effect, independently of the details of the flow and sediment transport models, which is 
by no means consistent with the formation of river dunes by a primary instability. In 
the last section, we present field experiments performed in the Leyre river, in France, 
that directly evidence the formation of ripples by a linear instability and the formation 
of dunes by a non- linear pattern coarsening limited by the free surface. Further obser- 
vations suggest that mega-dunes form when the sand bed presents hetereogeneities such 
as a wide distribution of grain sizes. Finally, we review and discuss the data available in 
the literature in the light of these theoretical and experimental results. 



2. The saturation length paradigm 

2.1. Saturation transient 

We present here a general framework to describe, with only few key variables, the sedi- 
ment transport properties. This framework can accomodate various situations associated 
to different modes of transport (bed load, saltation, reptation, etc) and different dynam- 
ical mechanisms controlling this transport (hydrodynamical erosion of the bed, splash, 
mixing by turbulent fluctuations, etc). 

The reference situation is a uniform turbulent boundary layer of constant shear velocity 
over a flat sand bed characterised by a threshold shear velocity Uth- In this situation, 
one observes a steady uniform transport characterised by a flux g sat called the saturated 
flux, which corresponds to an equilibrium state between flow and transport. g sat is a 
function of and Uth- It is important to emphasise the precise meaning of in this 
context. On the one hand, u* is defined from the Reynolds shear stress, i.e. from the 
averaged correlation of velocity fluctuations u'w'. On the other hand, the transport is 
mostly determined by the velocity field inside the transport layer of thickness ho (few 
grain diameters in the case of bed-load and typically 1 cm in the aeolian case) . If the 
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latter is embedded in the inner layer, then the average velocity is locally logarithmic 
and can be related to the basal shear stress, should thus be determined from velocity 
measurements performed inside this inner layer - i.e. typically at a distance of I ~ 10 -2 A 
from the bed, see part 1. The time-scale over which the momentum is exchanged increases 
linearly with the distance z from the surface (it scales as z/u*). Thus, the velocity 
fluctuations at time-scales smaller than z /u* should be associated to the turbulent shear 
stress and those slower than z/u*, to variations of u*(t). The shear velocity u*(t) relevant 
for sediment transport should thus be averaged over the time-scale ho/u*. 

Considering an inhomogencous situation - for instance a spatial variation of the shear 
stress - the sand flux is not instantaneously in equilibrium with the local shear stress. 
In most of the situations, the transient toward equilibrium can be described by a first 
order relaxation law, with a single time and length scales: 

T sat d t q + L sat d x q = q sat - q, (2.1) 

where the flow goes in the increasing x-direction. In a situation homogeneous in space, 
T sa t is the time needed for the flux to reach again the saturation if the flow speed suddenly 
changes. Conversely, in a steady situation where there is no sediment flux at the entrance 
of a volume of control, L sa t is the length needed for the flux to reach g sat . 

As the relaxation time is usually much smaller than the time-scale of evolution of the 
bedform, we are left with a description of the transport by three variables: the satu- 
rated flux (7 sa t, the threshold shear velocity u t h and the saturation length L sat . Note 
in particular that gravity effects are included into the transport threshold. This multi- 
scale approach allows to separate clearly the dynamical mechanisms that govern the 
emergence of bedforms from those governing sediment transport. In other words, sub- 
aqueous ripples and aeolian dunes are of the very same nature: a primary instability 
that is not sensitive at all to the mode of transport (e.g. bed-load vs saltation, see 
Her sen et al. 20021 IClaudin fc Andreotti 2006]) . This framework lies on two important 
assumptions. First, the depth of the inner layer in which the shear stress is vertically 
homogeneous should always remain larger than the depth of the transport layer - see 
part 1. Second, there is in general more than one relaxation mode and consequently a 
whole spectrum of relaxation times and lengths. The description by a first order relax- 
ation implies that one of these modes is significantly slower than (and thus dominant in 
front of) the others. Otherwise, two or more saturation lengths have to be taken into 
account in a higher order relaxation equation. 

An important such case is the suspended load, which must be described by a field 
(the sediment concentration) and not by a scalar (the integrated flux). Of course, a 
saturated flux g sat exists in that case (the only particularity is the explicit dependence of 
9sat on the flow depth) but, there is an infinite number of modes describing the behavior 
around saturation. We will thus restrict the discussion to the case of grains sufficiently 
heavy to prevent suspension. This is realized in practice for ripples and dunes in rivers, 
at moderate Froude numbers. It has been observed by |Guy et al. 19 66 that antidunes, 



which form in above T = 1, produce so much turbulent kinetic energy that suspended 
load could be dominant. This situation, probably highly non- linear, and which requires 
a specific treatment fsee lParker 19751 [Colombini 20041 IColombini fc S tocchino 2005}), is 
beyond the scope of the present article. 

We aim here to discuss the origin of the existence of a transport saturated state and 
to propose a self-consistent model, predicting the main characteristics of transport: the 
static threshold, the saturation length and the saturated flux. 
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2.2. Static threshold 

The bed load threshold is directly related to the fact that the grains are trapped in 
the potential wells created between neighbouring grains at the sand bed surface. To 
obtain the scaling law of the threshold shear stress, the simplest geometry to consider is 



a single spherical grain jammed between the two fixed grains below it (see figure 18 in the 
Appendix [A}. We consider the situation in which the cohesive forces between the grains 
are negligible and the friction at the contacts is sufficient to prevent sliding. The loss 
of equilibrium occurs for a value of the driving force F proportional to the submerged 
weight of the grain: F cx (p s — pf)gd 3 , where p s is the mass density of the sand grain, and 
pf is that of the fluid. As F is proportional to the shear force rd 2 exerted by the fluid, 
the non-dimensional parameter controlling the onset of motion is the Shields number 0, 
which characterises the ratio of the driving shear stress to the normal stress: 

e = ? PfU \ , (2.2) 

\Ps - Pf)ga 

The threshold value can be estimated from a force balance. In the viscous regime, 
the grain drag coefficient Cd is inversely proportional to the grain Reynolds number 
1Z = Ud/v. Cd = s 2 lZ~ 1 ; it is constant {Cd = Coo) in the turbulent regime. Note that 
the turbulent drag is not only due to the fluctuations induced by the grain itself, as in 
the case of a grain falling in a fluid at rest, but also to those present in the turbulent 
flow. This situation has never been studied properly so far. Note also that the drag 
force is a priori different in the vicinity of the ground and that the lift force could be 
non-negligible when the grain is trapped at the surface of the sand bed. In order to avoid 
introducing two many parameters, we will not consider these effects here. 

Typical values found in sedimentation experiments performed with natural sand grains 
are Coo — 1 and s ~ 5. The generic situation is a uniform shear stress p/u 2 , for which 
the velocity profile is linear in the viscous regime: 

u x = ^ (2.3) 

and logarithmic in the turbulent regime: 

^--ln(l + 4), (2.4) 
k \ rd) 

where k is the von Karman constant and r the ratio of the aero or hydrodynamic rough- 
ness to the grain diameter. The threshold Shields number is constant (see Appendix [A|, 
both in the viscous 



e,„ = 5 (2.5) 



and in the turbulent regimes: 



3Cooln (1 + l/2r) 

where p is the avalanche slope. The laminar-turbulent transition, indicated by the grey 
zone in figure |2j;, is determined by the Reynolds number and thus depends on the grain 
diameter d. In most of the situations relevant in geophysics, d lies in the cross-over 
between these two regimes. Figure [2] shows the comparison between experimental data 
for natural sand grains and the full model, whose derivation is provided in appendix [A} 
There is no adjustable parameter in the model as one can independently measure the 
avalanche slope p — tan 32°, the drag coefficients and the sand bed roughness r = 0.1. 
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Figure 2. (a) Threshold Shields number Q as a function of the grain size d, for nat- 
ural sand grains in water. Symbols: measurement by lYalin fe Karahan 19791 (o) and by 
Fern andez Luque fc van Beek 1976] (■). Solid line: Comparison with the model (equa- 
tions |A 7| and |A 8[ ). (b) Dependence of the threshold Shields number 0" h on the bed inclination 
a. Symbols: measurements by |Fernandez Luque fc van Beek 1976| for natural sand grains (■) 
and by Loiseleux et al. 2005 usi ng gl ass beads (o) in a narrow channel. Solid line: expectation 
of the model (equations A 7 and A 8 1, for natural sand grains. Dotted line: approximation by 
cos a + sina/fJ,, with /j, = tan 32 . (c) Threshold shear vel ocity Uth as a function o f the grain 
size d, for natural sand grains. Symbols: measurement by lYalin fc Karahan 19791 (o) and by 
Fernande z Luque fc van Beek 1976] (M). Solid line: Comparison with the model. Dotted line: 
asymptotic behaviour in the turbulent regime. Dashed line: asymptotic behaviour in the viscous 
regime. Grey zone: zone in which the viscous sub-layer is larger than the bed roughness: for 
grains smaller than d = 100 /im, the bed is hydraulically smooth for any shear stress u* above 
the threshold u t h- 



It is worth noting that continuum mechanics cannot predict correctly this threshold. 
Assuming that the granular presssure P and the shear stress r are continuous at the 
fluid/grain interface, the standard Coulomb criterion t/P = p leads to a null threshold 
shear stress. Part of the problem comes from the fact that the pressure P felt by the 
first layer of grains is due to its own weight: ~ (p s — Pf)gd. Introducing such a pressure 
discontinuity, one finds a threshold Shields number O t h — p, which is independent of 
the flow regime and much too large, compared to measurements. This means that there 
is no continuity between the fluid and the granular shear stresses, making continuum 
descriptions of sediment transport inherently problematic. 

2.3. Gravity or slope effect 

Due to gravity, it is much easier to transport sediments along the lee side of a sand 
bump than on the stoss slope. This effect is directly encoded into the dependence of the 
saturated flux on the transport threshold, which depends on the slope. Let us note a 
the bed angle at the scale of the saturation length. The force balance now reads: 

^(Ps - Pf)gd 3 (sine* + ^cosa) = |c d p/ (U& df . (2.7) 

and is the fluid velocity around the grain at the threshold (see Appendix |A| . We 
keep the notations w t h and 6 t h for the threshold shear stress and the threshold Shields 



9 



number for an horizontal bed (a — 0) and denote by u" h and 0" h the thresholds in the 
case of an inclined bed. 

The dependence of the threshold on the slope has been measured experimentally by 
Fernandez Luque fe van Beek 1976] for sand grains, in the turbulent regime. The model 



matches quantitatively the results, without any adjustable parameter (solid line on fig- 
ure |2Jd). Neglecting the fact that Cd is a function of the grain Reynolds number, the 
threshold shear stress can be written as: pj (u" h ) 2 = PfU 2 h( C0S a + sina//i). Inserting 
this expression into the saturated flux formula, one can see that gravity leads to a down- 
slope component of the flux, i.e. to a diffusion mechanism. This approximation still gives 
a reasonable fit of the data (dotted line in figure [2j)) . 

The slope dependence of the threshold shear stress has also been measured in the 
case of aeolian transport, for rough sand grains (Rasmus sen et al. 1996[) and is again 
perfectly fitted by the model. However, there exists an unexplained discrepancy in the 
case of spherical glass beads in a viscous liquid (Fig. J2Jd), as one would need a very low 
effective friction /x ~ tan 70° much lower than the avalanche slope (/x ~ tan 24°) to fit the 
data (fLoiseleu x et al. 20 05 ) . With such a value, the gravity effect becomes completely 
negligible. 

Importantly, the whole gravity effect on sand transport can be encoded into the thresh- 
old slope dependence, at the linear order in a. This was directly proved experimentally by 
Fernand ez Luque fc van Beek 1976] in the sub-aqueous case and bvfRasmuss en et al. 1 996 
in the aeolian case. In different articles (Richards 1980 and followers), /x was considered 
as a parameter that can be adjusted to fit the initial ripple wavelength. As they were 
missing the correct stabilising mechanism, aberrantly low values of the repose angle were 
found. In section 3.4 (see also Fig. [5]), where the wavelength selection issue is discussed, 
we shall keep for /x the avalanche slope determined experimentally. 



2.4. Erosion limited transport 

When a sand bed is submitted to a flow, only a small fraction of the grains at the 
surface are entrained and this erosion process takes some time to occur. We first model 
the erosion rate, and in particular the role played by the disorder of the surface. We 
define 7V(0) the fraction of grains at the surface susceptible to be entrained at a Shields 
number O. In the absence of flow, for — 0, M is null. It increases very quickly 
around the mean threshold Shields number t h and reaches 1 for large velocities. 7V(0) 
reflects and encodes the distribution of potential wells at the sand bed surface. Whatever 
the functional form chosen for this quantity, there are two important parameters: the 
Shields number around which J\f switches from to 1, which reflects the overall transport 
threshold, and the typical slope N' in the same zone. A simple choice is: 

A/"=0 if 0<0 m , 

N = rf ~ ^ if ©m<0<0«, (2-8) 
OM — o m 

AT = 1 if 0>0 M . 

For the sake of simplicity, we consider that one grain at the surface occupies an area d 2 . 
We start from a very simple assumption: when one grain is entrained, we hypothesise 
that there can be no further erosion in the surrounding area until the grain has moved 
by a distance comparable to its own diameter d. The description of the elementary jump 
by one grain diameter is derived in appendix [B] The erosion time T is the time needed 
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to move over a distance d. It behaves as: 



T oc In 



p f d 



in the turbulent regime and 



T oc In 



e + e 



th 



e-e 



th 



3/2 



(ps - Pf)g 



(2.9) 



(2-10) 



in the viscous regime. This means that the mean velocity goes to at the threshold like 
1/ In [0 — t h] , which corresponds to a bifurcation much sharper than the standard sad- 
dle node bifurcation in \/8 — 6 t h- The experimental observation (lAbbott &; Francis 19771 
Fer nandez Luque fc van Beek 1971)} ICharru et al. 2004|1 that the grain velocity - i.e. the 
mean velocity conditioned by the fact that the grain moves significantly at the time reso- 
lution of the instrument - is non zero just above the threshold is thus fully consistent with 
this expression. Far from threshold, the grain velocity is just proportional to the flow ve- 
locity (lAbbott & Francis 1977] |Fernandez Luque fc van Beek 19761 ICharru et al. 2004|) . 

The erosion flux <p^ is the volume of grains eroded per unit time and unit surface. It 
can be expressed as an integral over all possible local configurations: 



-de 



thi 



(2.11) 



where 



ird 3 f e AA'(9 th ) 
60 J Q d 2 T(e th ,6)' 

is the bed volume fraction. A good approximation of this integral is given by: 



tpi cx 



(e - e m ) 



4> (e M - e m ) 



111 1 e-e„ 




i)gd for e m <e<e M , (2.12) 



where c m ~ 1.41 in the turbulent regime and c m ~ 0.605 in the viscous regime. In 
particular one can recover from this expression the expected asymptotic behaviours. 
Close to the threshold, there is a regime in which the density of mobile grains is not 
sufficient to induce a significant reduction of the flow strength in the transport layer, and 
the grains can be considered as isolated, which corresponds to a dilute transport layer. 
The saturation of the flux is then controlled by the erosion rate and the disorder of the 
sand bed. This idea is due to ICharru et al. 2004[ who proposed a semi-phenomenological 
model that uses extensively experimental results. In the present model, the length of a 
trajectory is independent of the trap in which the grain was initially at rest and fluctuates 
around a well defined average value C The mean area explored by the grain is Cd and 
contains, by definition of C, a mean number of potential wells sufficiently deep to trap 
the grain equal to 1. We then obtain the expression: 

£ = l4(e)' < 213 » 

At small O, L is simply one grain diameter d; at large 0, Af tends to 1 so that C 
diverges. This divergence has been directly evidenced experimentally in the viscous case, 
by ICharru et al. 20081 In that case, the measured value of Qm is around 20 m . C is the 
length needed for the transport to adapt to a change of shear stress. It thus gives the 
saturation length L sat . 

This behaviour constrains the domain of validity of this regime: the shear velocity 
should be around the threshold, the transport should be intermittent and diluted and the 
distribution of potential wells should be very large. As suggested by the experiments of 



11 



20- 


Qsat 


(a) 

□ 


□ / 


10 1 " 
10°- 


15- 




□ i 


OA 

A D 


io-l 


10- 

5- 




a a Q A 

W° 

crnr a 
nOg □ 
MiD □ 


□ 

□ 
□ 


10- 2 _ 
10" 3 - 


0- 








io- 4 



0.2 0.4 0.6 



0.8 1 






10" 



Figure 3. Saturated flux g sat rescaled by \ j ( 







10 u 




10"' 



10-2 ^lO" 1 

0-0 th 



L) gd* 



as a function of the Shields number 



0: raw data (□) gathered by Julicn 1998; same after local averaging (•); raw data (o) obtained 
by Fernandez L uquc & v an Be ck 1976 The solid line is the best fit by the model described 
here (equations 2.17| and |2.18[ ), which gives: © t h = 0.045 and 9m = 5 6a. (a) Lin-lin plot; (b) 
Log-log plot; (c) Log-log plot as a function of — t h- 



ICharru et al. 20041 a granular bed prepared by sedimentation is initially very disordered 
and consequently presents a wide range of potential wells. It then takes a very long time 
for the surface to re-arrange, leading to a drift of the distribution Af toward larger and 
larger threshold shear stresses. This long transient is probably not relevant to geophysical 
situations in which the time and length scales of the systems are always sufficiently large 
to ensure that an equilibrium state for the geometrical arrangement of surface grains has 
been reached. The saturated flux is simply the erosion rate times the hop length C: 



9sat oc 



(e - e m ) 



<t> (0m - e) 



•In 



e+e„ 
e-e„ 




-i)gep for e m < e < e M . (2.14) 



By construction, it vanishes at the Shields number m and diverges at @m- The trans- 
port is thus limited by the erosion rate only close to the threshold O m , in the presence of 
a disordered bed. As the density of transported grains becomes important, the transport 
becomes limited by the available momentum. 



2.5. From erosion limited to momentum limited transport 

One can extend the previous model by including the negative feedback of the transport 
on the flow ( TBagnold 19 56). Each time the flow entrains a grain from the bed and 
accelerates it, the grain exerts in return a stress on the fluid. The balance between 



erosion and deposition is the same as in the previous paragraph (Eq. 2.141 except that 
the shear stress in the transport layer is reduced to some value Uf, compared to that in 
a particle-free flow: 



<?sat oc 



"thj 



<t> Om - u \ ) 



In 




1 gd\ 



(2.15) 



We consider the regime in which the grains are transported at the surface of the bed i.e. 
do not form a surface sheet flow. They leave the bed with a velocity v-\ and collide back 
the sand bed with a velocity vi . The sand-borne shear stress is proportional to the sand 
flux and to the difference (vj — «t). The fluid in the transport layer is assumed to be 
in equilibrium between the driving shear stress PfU^, the fluid-borne basal shear stress 
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Pfuf and the sand-borne shear stress: 



Pf u l = Pf u \ + Psfi 



("1 



r) 



c 



(2.16) 



This relation is nonetheless valid in the saturated state (q = (feat) but also during the 
transient of saturation. It is very important to understand that the particles are entrained 
by the flow at a velocity reduced by the presence of other particles. The grain trajectory 
is encoded into a single quantity, (v± — v^)/C, which is a function of the reduced shear 
velocity u{ and not of it*. For the sake of simplicity, we limit the discussion to the 
turbulent case. The transposition to the viscous case is straightforward. Assuming that 
U| — scales on the grain mean velocity, the ratio (v± —v^)/C is proportional to the hop 
time T, which is also the erosion time. Equation (2.161 can be solved to get the flux: 



Pi i i 2\ rr/ x Pf Pf d 

9sat OC — («;-Uf) T{Uf) oc — a i — 

Ps p s V yps - Pf)g 



111 



U{ + u t h 



u { - u t h 



(2.17) 



Eliminating the flux between the equations (2.171 and (2.151, one obtains a relation 
between itf and u*. 



^ (^-l)gd 



4> (u^ - uf) 



111 



In 



Mf+«th 
Uf-U th 



Just above the threshold, the flow in the transport layer is undisturbed: Uf ~ u 
erosion limited regime is recovered and the flux can be approximated by: 



9sat oc 
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1 gd* 



(2.18) 
t . The 

(2.19) 



Far above the threshold, for large values of it* , the shear velocity it / inside the transport 
layer tends to um i-e. to a value independent of it*. The flux then scales as: 



9sat oc 



p f d 



PJ_ 

Ps V (p s ~ Pf)g 



(2.20) 



Thus, in the momentum limited regime, q scales at large velocities as u\ and not u\ as 
usually obtained. The same result is valid for aeolian transport, for which the negative 
feedback of particles on the flow has been directly evidenced experimentally. Although 
many authors have followed Bagnolds, the scaling law in ul is, in that case, the best 

we have re-plotted the mea- 
for the subaqueous case. It 



description of existing data (lAndreotti 2004)1 . In figure [3 



surements of the saturated flux collected by Julicn 1998 



exhibits a large dispersion due to several factors. First, data coming from systems with 
different grain size distributions have been plotted together without any distinction of 
symbols. Second, the reproducibility of experiments is made difficult by the problem of 
granular bed preparation: we emphasise again that sedimentation leads to an out of equi- 
librium situation that can last for days. Third, it is difficult to estimate the basal shear 
stress in flume experiments due to lateral boundaries. We have also plotted the measure- 
ments performed by |Fernandez Luque fc van Beek 1976[ which are much less scatterred. 
Given these reservations, the fit by the model derived here provides a good description 
of existing data. In particular, the asymptotic behaviors close to the threshold and for 
large shear velocities are well captured. 

As many other formulas would provide an equally good fit (i.e. |Meyer-Peter fc Miiller 1948] 
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Ein stein 19501 |Bagnold 1956| or lYalin 1963| . forthcoming sediment tranport dedicated 
experiments should be designed to test directly the ingredients of the models. For in- 
stance, the measurement of the hydrodynamic roughness in presence of transport can pro- 
vide a direct evidence of the negative feedback of particles on the flow fsee lAndreotti 2 004 
for the aeolian case) . One expects this roughness to be larger than in the transport-free 
case as a consequence of the increase of the dissipation rate - the fact that the particles 
are moving does not lead to a drag reduction. 

The last key issue is the saturation length. The sand flux is the product of the grain 
velocity by the density of mobilised grains. Two important mechanisms can thus control 
the saturation length: the length needed by the grain to reach its asymptotic velocity 
- the so-called drag length - and the length needed for erosion to take place i.e. the 
trajectory length C. The modelling of the drag length is a difficult problem as the 
trajectory takes place in a turbulent flow whose fluctuations are not due to the motion of 
the grain itself. The problem is thus very different from that of a sphere moving in a fluid 
at rest, a problem for which the drag law is calibrated. To the best of our knowledge, the 
motion of a sphere whose diameter lies in the inertial range of the turbulent flow is still 
an open problem. We are thus left with the standard drag force formula ^CdPfU 2 d 2 , 
with a drag coefficient of order one. Then, solving the equation of motion, one obtains 
a drag length Ldrag around 2 (p s /pf)d. In summary, we consider here that L sat can be 
either limited by erosion, in which case it is expected to scale as: 

iaat - d (2.21) 

or by grain inertia, in which case it should scale as: 

L sat ~ 2^d (2.22) 
Pf 

One can see that these two predictions are difficult to test and discriminate. The ero- 
sion length gently increases with the shear velocity while the drag length is independent 
of it; the drag length increases with the density ratio p s /pf but, experimentally, this 
parameter cannot be easily varied by a large factor. We shall see below that the study 
of subaqueous ripples can shed some light on this issue. 
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3. Subaqueous ripples 

3.1. Linear stability analysis 

Given that the details of the sand transport only affect the saturated flux, the saturation 
length and the threshold shear stress, we can perform the linear stability analysis of a 
flat sand bed in a very general framework. We consider a periodic disturbance of the 
bed profile Z. As the base state is homogeneous, we can seek for modes of the form 
exp(at + ik(x — ct)). The shear stress a xz induced by the wavy sand bed is written in 
the Fourier space as 

a xz = ul(A + iB)kZ, (3.1) 

where it* is the shear velocity in the reference state and k = 2t:/X the wave number 
associated to the spatial coordinate x. We refer the reader to the part 1 of this article for 
the derivation of the components of the shear stress A and B respectively in phase and 
out of phase with respect to the bottom profile. The lag between the maximum shear 
stress and the ripple crest is given by \B/(2irA). The computation of the coefficients A 
and B has been performed under the assumption that the sand bed can be considered 
as static. As a matter of fact, the growth rate a is related to the sediment transport 
and on the order of g sa t/i| at . Experimentally, a is usually found to be four orders 
of magnitude smaller than the typical flow shear rate: the time-scale of formation of 
subaqueous ripples is typically ten seconds to be compared to few tens of microseconds 
for the period at which the flow is excited by the dune relief, X/U. Therefore, in contrast 
to what is suggested by IColombini fc Stocchino 20051 the normal velocity of the grain- 
fluid interface does significantly change the flow in this problem. 
The threshold shear stress is a function of the slope: 

"th = w th-^-, (3.2) 

where (i is the avalanche slope. At linear order, we can write the saturated flux q sat (u* , u t h) 
as: 

<7sat = -5 —{A + iB)kZ+- —ikZ (3.3) 

cm* 2 cm t h 2/x 

Introducing a reference flux Q, we rewrite this expression under the form: 

?sat = Q(a + ib)kZ (3.4) 

where a and b are the components of the saturated flux in and out of phase with the 
topography. Now the flux relaxes to its saturated value with a spatial lag L sat : 

ikL sat q = q sat - q (3.5) 

Using the conservation of matter d t Z + d x q = 0, one obtains the dispersion relation: 

., ., Q ik g S at iQ(a + ib)k 2 

a — ike = —ik— = —- — — = — (3.6) 

Z 1 + ikL sat Z 1 + ikL sat 

Splitting the equation into its real and imaginary parts, one obtains the growth rate a 
and the propagation speed c: 

■^satc = (kL sat ) 2 (b- akL sat ) 

Q 1 + (fci sat )2 1 • ] 

L sat c _ (fcl/ sat )(a + bkL sat ) 

Q ~ 1 + (fcisat) 2 1 ' ' 



15 



This corresponds to a standard convective instability at large wavelengths (Fig.|4|. The 
cut-off wave-number k c above which modes are stabilised by the saturation length is 
given by: 

k c L sat = - ■ (3.9) 
a 

Note that the instability can present a different threshold that that for the transport if 
b vanishes at some value of larger than u t ^. In first approximation, a and b are weak 
functions of kz . Then, one can approximate the maximum growth rate wave-number 

km ax aS . 



x -l/3_ x l/3 with x 




(3.10) 



Its 



As the instability is convective, we have also computed the spatial growth rate 
maximum nicely coincides with that of the time growth rate. 

It can be inferred from the conservation of matter that the propagation velocity c is 
proportional to the difference of flux Sq between trough and crest and invertly propor- 
tional to the ripple height 2£ - this is the so-called Bagnold relation. At large wavelength 
A, 5q is proportional to the reference flux Q and to the height so that the propagation 
speed varies as: c oc Q/X. This is confirmed by figure |4Jd, which shows a roughly linear 
relation between c and k. 

3.2. Erosion limited transport 

We can make the previous general results more precise by making explicit the expressions 
for the length L sat and time L 2 at /<3 scales, and for a articular relation between sediment 
flux and shear stress. In the case of the erosion limited transport, the saturated flux 
reads: 



<?sat = X- 



- In 1 



(3.11) 



Defining the reference flux as 



Q = x 



In 1 



we get: 



a = A 



1 



b=B 



(uth/u*) 2 

ln(l-( Uth /«*) 2 ) 

(uth/u*) 2 

ln(l - (uth/w*) 2 ) 



1 



(3.12) 



(3.13) 



1 



ln(l-( Wth /^) 2 ) 



(u th /u*) 2 
A* 

Close to the threshold, (for it* ~ itth), a tends to A and b to B — . For asymptotically 
large shear stresses i.e. ^S> u t h, a tends to A and b to B. Thus, the ratio b/a increases 
with the flow strength. This means that the destabilising mechanism becomes more 
efficient as increases so that the maximum growth rate wavelength A max decreases 
with (Fig. [6). 

If B is below then b vanishes for a value of it* larger than uth- The instability 
threshold is then distinct from the transport threshold, as observed in the viscous regime. 
As expected for such an instability, the most unstable wavelength then diverges at the 
instability threshold (Fig. [6]d) - in practice, it would be limited by the geometrical size 
of the experiment. In the turbulent regime, however, B is sufficiently large to ensure that 
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Figure 5. Wavelength A max of maximum growth rate as a function of the ratio of the two 
characteristic lengthscales zo/L sa , t , in the limit u, 3> ti t h, for [i = tan 32°. Black symbols: ratio 
A ma x/isat- Open symbols: ratio A max /zo- Squares: erosion limited model. Circles: momentum 
limited model (One can hardly see differences between the two transport models). 



flat beds are unstable as soon as transport takes place. Correspondingly, A max remains 
finite at = u t h- 

3.3. Momentum limited transport 

When the transport is limited by momentum, we can write the saturated flux under the 
generic form: 

<7sat=rf(^-u t 2 h ). (3.14) 



The self-consistent model derived above gives 7 = 0, while the Bagnold 1956 formula in- 



volves an exponent 7 = 1/2. Other empirical models such as |Meyer-Peter fc Miiller T 948 , 
Einstein 1950 or Yalin 1963 can also be approximated in this way. We now define the 
reference flux Q by: 

Q = ( 7 + 1)X^ (7+1) , (3.15) 



so that a and b read: 



A - 



1 A < 
1 + 7 
jB + fj,- 



b = B-^— t f. (3.16) 

1+7 

For asymptotically large shear stresses i.e. u* ^> u t h, a tends to A and b to B, indepen- 
dently of 7 and fx. Close to the threshold, (for u* ~ u t h), a tends to A/(l + 7) and b to 
(B — /i _1 )/(l + 7). So, again, the instability threshold is above the transport threshold 
when B is smaller than [iT 1 . 

3.4. Ripples wavelength selection 

Using the results of part 1, the fastest growing wavelength A max can be computed, taking 
into account the dependencies of A and B on kzo. The most important question is the 
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Figure 6. Wavelength A max rescaled by L sa t as a function of the rescaled shear velocity M*/«th, 
for jj, = tan 32° and L sa t/zo = 80 (typical subaqueous case with sand grains). Erosion limited 
model (■), momentum limited model (7 = 0, o) and Bagnold 1956] model (7 = 1/2, •). (b) 
Same plot, computed with the erosion limited model for different values of [i: tan 24° (o), tan 32° 
(■), tan 70° (□). 



scaling law followed by this wavelength. There arc two length-scales in the problem: z 
which is an hydrodynamical quantity and L sat which is related to the sand transport. 
Figure [5] shows that the saturation length controls the scaling of A max if Asat is larger 
than zq, whatever the transport model. Conversely, zq controls the scaling of A max if it is 
larger than 10 L sa t- This could be the case in the hydraulically smooth regime for which 
zq scales on v/u* QSumer fc~ Bakioglu 1984[ ). For subaqueous ripples, the measurements 



of initial wavelength are usually larger than 100 d (see Fig. 10), while Zq is on the order 
of 0.1 d. If the scaling was controlled by z$, one would underestimate the most unstable 
wavelength by two orders of magnitude. As a consequence, the wavelength is controlled 
by the transport saturation length L sa t. This means that models in which a univoque 
relation between the actual flux and the shear stress is used cannot capture correctly 
the physics of ripple instability. The second conclusion is that, although the saturation 
length may be determined by different dynamical mechanisms, aeolian dunes are of the 
same nature as subaqueous ripples. It means that different modes of sediment transport 
(saltation, reptation bed-load, etc) in different situations (viscous, turbulent, etc) can 
lead to bedforms instabilities of same nature]]] 

Figure [6^ shows that the transport model as a negligible influence on the selected 
wavelength. This evidences that L sat and q S at are the single two relevant quantities 
encoding the sediment transport details. Due to the stabilising role of gravity encoded 
in the slope dependence of the threshold, A max increases close to the threshold shear 
stress. Of course, this effect is very sensitive to the value of fi, as shown in figure [6)3. 
For natural sand grains (/i = tan 32°), the wavelength A max decreases from 30 L sat at 
the threshold to 20 L sat far from it. For glass beads, the prediction is very different 
depending whether one takes the avalanche slope (/1 = tan 24°) or the experimental data 
of ILoiseleux et al. 2 005 (equivalent to a fictitious /1 = tan 70°) for the slope effect. In 
the later case, A max is almost independent of m* whereas it diverges at the threshold in 
the former case as B is then of the order of l//i. 

f Of course, the aeolian ripples do not belong to the same class of bedforms as they result from 
a screening instability ( |Bagnold 1941] lAnderson 19871 lAnderson 19901 lAndreotti fc al. 2006). 
not from a hydrodynamical instability. 




Figure 7. Ripples aspect ratio 2£nl/A selected by hydrodynamical non-linear effects as a 
function of the wave- number k rescaled by L sat . The dotted line corresponds to a calculation 
based on the momentum limited transport model with L sat /zo = 350 (typical aeolian case), 
u*/uth = 3.8, /i = tan 32° and 7 = 0. The dashed line corresponds to the erosion limited 
transport model with the same parameters except L sa ,t/zo = 80 (typical subaqueous case). The 
solid line does not take any saturation length into account and serves as a reference: the sharp 
decrease of the two other curves are associated with the stabilisation of small wavelength by 
the saturation len gth. The symbo ls correspond to subaqueous measurements performed by 
Baas 1999] (•) andfGuy et al. 1966|(A) in flume experi ments and to aeolian field measurements 
performed bv lParteli et al. 20061 IBaddock et al. 20071 (D). 



3.5. Ripple amplitude selection 

After the linear stage during which the ripples emerge, they exhibit pattern coarsening 
by progressive merging of bedforms. During this slow process, each of these ripples may 
be considered as purely propagative structures that do not grow nor shrink. In the course 
of their evolution, they should thus present an amplitude 2£ which is only function of 
their wavelength A. 

In part 1, we have performed a weakly non- linear expansion of the hydrodynamics 
above a topography up to the third order in kQ, which is the small parameter of the 
problem. As a main output of these non-linear calculations, the phase shift between the 
elevation profile and the basal shear stress decreases as the amplitude Q is larger, which 
stabilises the bedform. Considering a sinusoidal profile, the Taylor expansion of the basal 
shear stress components takes the form A = A\ + A^kC,) 2 and B = B\ + B^(kQ 2 . We 
denote Cnl, the particular amplitude for which a purely propagative solution is selected, 
i.e. which does not grow nor decay: dtZ = —cd x Z. This gives the condition: b = akL sat - 

Figure [7] shows the prediction of the bedform aspect ratio 2£nl/A as a function of 
kL sat for two sets of parameters, which corresponds to typical aeolian and subaqueous 
situations. At small wave-number, the aspect ratio increases logarithmically with k. It 
then decreases abruptly and vanishes at the cut-off wave- number k c . For the aeolian 
conditions, the predicted dune aspect ratio is around 1/15 which quantitatively matches 
field measurements (~ 1/13). For current ripples, the measured aspect ratio are much 
more dispersed: Allen 1985 and Coleman & Melville 1994 report values around 1/20 
whereas [Baas 19 99 found larger values around 1/8. These discrepancies may be due to 
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Figure 8. Saturation length in water determined from experimental measurements of the initial 
wavelength as a function of the shear velocity for different types of particles, (a) Glass beads 
(Valance & Langlois 2005): d = 250 jj,m (•) and d — 500 /im (□). (b) Natural sand grains 
I Coleman fc Melville 19961 and iBaas 1999(1 : d = 210 f im (• ), 240 ^m (♦), d = 830 (□). 
The solid lines correspond to the best fit of equation (2.13l. These curves tend to diverge for 
U* = it a/ — 4.5 Mth or equivalently for G = Om ~ 20 Bth- In (a) and (b), the small factor 
between small and large grains could be due to a subdominant dependence of L sa t on viscosity. 



the measurement techniques related to the small amplitude of subaqueous bedforms (few 
tens of grain diameters) . 

3.6. Comparison with experiments 

Rather than predicting the wavelength and comparing it with experimental data, we can 
invert the process and determine the value of the saturation length that would give the ob- 
served wavelength. To play this game, we have chosen the erosion limited transport model 
and checked that other choices do not change the conclusions reached here, fi has be cho- 
sen equal to the avalanche slope. Most of the experimental data available in the literature 
correspond to well developed ripples. As these bedforms exhibit pattern coarsening, it is 
very important to focus on papers reporting the initial wavelength (linear regime), mea- 
sured for grain sizes larger than 200 /jm (hydaulically rough sand bed) . We have selected 
five such data sets: Colema n fc Melvi lle 1996, Baas 1999, Vala nce" fc Langlois 2005[ 

As shown in figure [8] they present consistent trends, which is not the case for the data 
obtained with smaller grains (d ~ 100 ^m). The saturation length is found to be on the 
order of several grain diameters. It is slightly smaller for the glass bead experiment than 
for the natural sand grain ones. The data points around the threshold are very sensitive 
to the values taken for /x and Uth in the model and should not be over-interpreted. The 
slight increase of L sat with is more robust, although it is based on the last few data 
points. More significant is the decrease of the ratio L sat /d when d increases, a feature 
present for both glass beads and sand grains. It could be related to a subdominant 
dependence on viscosity. 

In the section devoted to the description of sand transport, we have discussed the two 
simplest possibilities for the dynamical mechanisms limiting the transport saturation: 
erosion and inertia. Concerning erosion, the prediction L sat = C oc d/(l — Af(Q)) is 
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Figure 9. Comparison between the saturation length, rescaled by the drag length, in the 
aeolian (*) and subaqueous (other symbols, see Fig. [8| cases. For aeolian dunes, L sa t is deter- 
mined from the most unstable wavelength under well characterised winds Andreotti ct al. 2008 
For subaqueous ripples, each point of this graph corresponds to the very same data as in fig- 
ure [8] but averaged over 6 measurements. The dashed line corresponds to the entrainment 
threshold. The solid line is the average over the different points measured in the aeolian case: 
L sat ~ f .66 {ps/pf) d. 



almost impossible to verify as the distribution of potential wells at the surface of the 
bed is not known. Still, with the simple parametrisation chosen above, one expects the 
saturation length to scale on the grain size and to gently increase with it*. The solid 
line in figure [8]d is the best fit by such a form. One obtains the estimate: OM/©th ~ 20. 
This is far above the value found in figure [3j around 5, and that found experimentally 
by ICharru et al. 20081 around 2. It means that the erosion length is probably not the 
mechanism limiting saturation. Otherwise, one would expect a much more rapid increase 
of the saturation length. Instead, one can observe that it is almost constant, with a 
subdominant, slow increase with it*. 

Figure [9] aims to test the second simple possibility: an inertia limited saturation 
length. In this case, one expects a scaling law of the form L sat cx (p s /pf)d sim- 
ilar to that observed for aeolian sand transport (Hcrscn et al. 2002, Andreotti 2004. 
IClaudin fc Andreotti 20061 lAndreotti k. Claudin 2007|) . We have used the measurements 
of the wavelength at which aeolian dunes form reported in Andreotti ct al. 2008 , obtained 
either in the field or using aerial photographs. Using the inversion method proposed here, 
the saturation length L sat has been obtained for different values of the wind shear velocity 
(Fig. [9]) and is mostly independent of it. Once rescaled by p s /pj d, L sat is of the same or- 
der of magnitude for both aeolian dunes and subaqueous ripples. In first approximation, 
all series of data, considered separately, are independent of u*. A better agreement with 
aeolian data is observed for glass beads and for large grains. A discrepancy by a factor 
of two is observed for small natural sand grains. As a possible interpretation, the grains 
roll on the ground during their phase of acceleration, which may lead to underestimate 
the length needed to reach the fluid velocity. One expects rough sand grains to be more 
sensitive to this effect. 

For a given type of grain, the extra-dependence of the initial wavelength on the grain 
diameter may be attributed to viscous effects, either on the saturation length or on 
hydrodynamics. To test this, we have plotted in figure [10] the measurements of the 
initial wavelength together with different hydrodynamical calculations, assuming that 
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Figure 10. Initial wavele ngth A in water as a funct ion of the shear velocity for different types 
of particles: glass b eads (|Valance fc Lan glois 20 05]): d = 250 um (•) and d — 500 /im (■); 
natural sand grains (Coleman & Melville 1996 and Baas 1999): d — 210 fim (•), 240 fim (♦), 
d = 830 fim (□). The lines show the prediction of the model, assuming that the saturation 
length is given by the drag length, with an hydrodynamically rough sand bed (solid line), with 
a viscous surface layer (dotted line, IZt = 125), and in an intermediate regime (dashed line, 
Tit = 10). 



the saturation length is limited by inertia: L sat is kept fixed and equal to 2 (p s / pf) d. 
As the shift between the shear stress and the topography gets reduced for a viscous 
surface layer (see part 1), the predicted wavelength is larger by more than one order of 
magnitude. Moreover, the coefficient B becomes smaller than 1/ fi so that the instability 
threshold, at which the initial wavelength diverges, is above the entrainment threshold. 
This applies to grain sizes much smaller than 100 /im (or to more viscous liquids), in the 
hydrodynamically smooth regime. Still, it should give the trend for the data analysed 
here: as the grains size decreases, the influence of viscosity increases and so does the initial 
wavelength. Note that this trend is in the same direction as the empirical observation of 
IColeman et al. 20031 that X/d cx cT - 25 



4. Effect of the free surface 

4.1. Dispersion relation 

In the previous section, both hydrodynamic and erosion aspects have been gathered to 
study bedforms (ripples) under an infinite depth assumption. However, rivers have a free 
surface and the water depth H is finite. We expect these two additional ingredients to 
significantly change the shape of the dispersion relation. 

Our goal is to relate the formation of ripples and dunes to the two relevant length scales 
H and L sat . In part 1, we have precisely computed the basal shear stress coefficients 
(A and B) in the case of a turbulent flow over a wavy bottom. The large wavenumbers 
are insensitive to the free surface (Fig. 14 in part 1). By contrast, A and B display a 
resonance peak around kH ~ l/J- 2 and have a divergent behavior as k — > 0. As a and c 



are directly related to A and B through equations (3.7 1 and (3.8 1, they both present the 
same features. In figure [TTfc , the growth rate is represented as a function of kzo for three 
different values of H/z . All curves collapse in the large-A: region on the dispersion relation 
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Figure 11. Growth rate a and propagation velocity c as functions of the wave-number k 
for different values of H/zo (H/zq = 5000, 10000, 20000) and the typical values of parameter 
T — 0.8, n — tan 32°, u t h/u, = 0.8 and L Bat /z = 80. In (a), a is rescaled by L sat and k by zq, 
whereas in (b), (c) and (d) all lengths are rescaled by the relevant length-scale H, so that all 
curves collapse in region where kH ~ 1. The dashed line represents the reference unbounded 
case presented in the previous section. In panel (c), due to the log-scale, the absolute value of 
a is displayed and the grey areas encode for negative values of the growth rate. 



computed in the reference unbounded case. In particular they exhibit a maximum for 
the same wave-number, which corresponds to the initial ripple wavelength. This means 
that the presence of the free surface does not influence the formation of ripples as long as 
H and L sat are well-separated length-scales. In a zone around kH ~ l/T 2 (see Fig.flTjj), 
the function a(k) presents, in comparaison to the reference case, a sharp dip which can 
be attributed to a resonance of gravity surface waves, independently of the transport 
issue (see part 1). As shown in figure 12 1, the width and the amplitude of this dip is 
very sensitive to the value of the Froude number. For small J 7 , the effect of the free 
surface is marginal and the dip is very small. As the Froude number increases, the dip 
becomes more pronounced so that the growth rate a becomes negative in an enlarging 
range of wave-numbers: the free surface stabilises wavelengths commensurable with the 
flow depth. Last, the semi-logarithmic plot of figure [TTJ; reveals the behavior of the 
growth rate in the small-kH limit: er(fc) tends to from below. This indicates that the 
very large wavelengths are also stabilised. In the intermediate range of wavelengths, a 
slight increase of the growth rate is observed. 
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Figure 12. Dispersion relation obtained for a ratio H/zq — 10000, fi = tan 32°, iitt/w* = 0.8 
and L aa x/zo = 80. (a) Growth rate a as a function of the wave-number for different values of 
the Froude number, (b) Propagation speed c. In both panels, the Froude number T is varied 
from 0.5 to 1.2 by increment of 0.1. 



The other output of the linear stability analysis is the propagation velocity of the 
pattern c(k) (Fig. Ill), which is also very sensitive to the resonance at kH ~ 1/T 2 . It 



presents a sharp maximum on the left of the resonance followed by a dip on the right of 
it. At this dip, c can become negative for sufficiently large Froude numbers {T > 0.7). 

4.2. Dune formation 

Gathering the different dispersion relations in the (J 7 , kH) space, one produces the stabil- 



ity diagram (Fig. 13). The central region delimited by the two marginal stability curves 
corresponds to the zone of unstable wavelengths (a > 0). Large k are stabilised by the 
saturation length, which explains the lack of dependence on the Froude number in this 
zone: both the marginal stability and the maximum growth rate curves are vertical lines 
in the diagram. This simply reflects the fact that ripples do not feel the free surface: 
they disturb the flow over a thickness of the order of the wavelength A i.e. much smaller 
than H . As already mentioned, the most unstable mode always corresponds to ripples 



(dotted vertical line in the large- k zone of figure 13 1. 

A second zone of stable modes is located at much smaller kH and is associated to the 
influence of the free surface. The bedforms excite gravity surface waves that propagate at 
a speed y/(g/k) tanh(kH) in the frame of reference moving with the surface particles at 
the velocity M SU rface- When the upstream wave velocity balances the downstream surface 
velocity, for 



the disturbances induced by the bottom topography accumulate. In these resonant con- 
ditions, the wave amplitude reaches a maximum (Fig. 20b in part 1). The waves are 
in phase for wave-numbers above the resonance and in antiphase below it (Fig. 20c in 
part 1). At the resonance, the free surface is in quadrature with the bottom, which 
tends to move downstream the point of maximum shear stress i.e. to stabilise the bed- 
form. Of course, as the influence of the free surface on the flow is localised over a typical 
distance A, its effect is more important as kH becomes smaller. In summary, there 
are two conditions under which the free surface can overcome the inertial destabilising 
effect: (i) around the resonance, since the standing wave amplitude is very large and 
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Figure 13. Stability diagram parametrised by the Froude number T and the rescaled 
wave-number kH. The marginal stability curves (a = 0), shown in solid lines, separates the 
unstable zone (white) from the stable ones (grey). The overall maximum growth rate (bold 
dotted line) is always reached for ripples. The local maximum of the growth rate resulting 
from the resonance of standing surface waves is shown in thin dotted line. The dashed line 
(c = 0) separates upstream (c < 0, dark grey) from downstream propagating bedforms. The 
other parameters are set to: H/zq = 10000, n = tan 32°, itth/w* = 0.8 and Lsat/zo = 80. 



(ii) in the limit of small wavenumbers, as H becomes much smaller than wavelength A. 
This is precisely what can be observed in figure |13| the sharp stable zone surrounding 
the resonant curve. For obvious reasons, this new minimum of the growth rate associ- 
ated to the surface wave resonance comes with a local maximum of a (thin dotted line 
along the restabilised zone). In the work of Richard s 19801 the latter has been asso- 
ciated to the formation of dunes by a linear instability. In many other linear analysis 
([Kennedy 19631 I Reynolds 1965} [Engelund 1970 [ |Freds0e 19 741 IColema n fc Fenton 2000} 
Colombini 2004, Colombini & Stocchino 2005), the prediction of ripples is missed for the 
reasons detailed above (the sediment transport relaxation is not correctly taken into 
account), i.e. the corresponding maximum in the dispersion relation is absent. As a 
consequence, they are left with a unique peak in the region of kH around unity, which is 
found here to be a secondary maximum. As explained below, we fundamentally disagree 
with the conclusion reached in all these papers that subaqueous dunes result from the 
linear instability of a flat bed. 

Let us briefly recall the basic reasons for which one usually associates the appearance 
of a pattern to a maximum of the growth rate. One considers that the initial condition 
Z(t = 0) is essentially flat, with some wide-band noise. Its Fourier transform Z(k, t = 0) 
then contains some energy in a wide range of wavenumbers k. In the linear regime, the 
surface profile reads: 

Z{k,t) = Z{k,t = 0)e CT(fc)t (4.2) 

If the distribution of initial amplitude is initially sufficiently flat and if the growth rate 
cr(k) presents a sharp absolute maximum in fc max , then a pattern dominated by the 
corresponding wavelength A max emerges, as this mode grows the fastest. In the present 
case, the amplitude of this secondary maximum close to the resonance is almost the same 
as the value of a at the same wavenumber in the unbounded case. Moreover, most of 
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the modes between the resonance and the ripple peak are in fact much more unstable: 
if a linear instability could be invoked, the amplitude of all these intermediate modes 
would eventually be larger than that of this local maximum. Furthermore, the ratio of 
the primary and the secondary maxima of the growth rate is on the order of (H / L sa ±) 2 
which is a large number: it is typically on the order of 10 4 for flume or small river 
experiments (d ~ 400 /im and H ~ 40 cm) and 10 6 for a large river (d ~ 400 and 
H ~ 4 m). For example, taking half a minute for the characteristic ripple apparition time 
(see Fig. 14 1), it would give ~ 3 days for the dune linear growth time-scale, i.e. much 
too large in comparaison to observations (see Fig. 151. For these reasons, this secondary 
peak in the dispersion relation cannot be associated to dunes. 

Let us contrast this subaqueous situation to that of aeolian ripples superimposed on 
aeolian dunes. As already mentioned, the instability mechanism of these dunes is of 
the same hydrodynamics nature as that of subaqueous ripples, i.e. comes from the up- 
wind shift of the basal shear stress with respect to relief. Aeolian ripples, however, are 
generated by a screening instability: the upwind face of ripples receives more impacts of 
saltating grains than the downwind face (Bagnold 1941 lAnderson 19871 fAnderson 19 90). 
Rapidely, non-linearities makes the ripple pattern saturate to a wavelength much smaller 
than that at which dunes emerge, and this saturation is faster than the dune time for- 
mation. One can then consider that dunes results from the linear instability of a flat bed 
which presents saturated aeolian ripples. Moreover, the growth rate at the wavelength 
of ripples in the dune instability is negative, because ripples are smaller than the aeo- 
lian saturation length. Conversely, the growth rate at the wavelength of dunes in the 
ripple instability is much smaller than the growth rate at the same wavelength in the 
dune instability. In conclusion, aeolian ripples and aeolian dunes can truly be associ- 
ated to two different linear instabilities. In the case of subaqueous ripples and dunes, 
none of these criterions (different destabilizing mechanisms, saturation of ripples wave- 
length, separation of the modes by several decades of almost non-growing wavelengths) 
is full-filled. 



5. Field experiments 

In this last section, we present direct experimental evidences that river dunes do not 
form by a linear instability. We will discuss our field measurements in the light of the 
model proposed here and show reciprocally that all existing observations are consistent 
with this model. We will finally propose a new definition of the different subaqueous 
bedforms, based on the physical mechanisms which control their formation. 

5.1. Formation of ripples 

We have studied the formation of ripples, dunes and mega-dunes in the Leyre river. This 
river is located in the south west of France, in a region called 'les Landes de Gascognc' 
(44°32'N, 0°52'W). It flows in a particularly homogeneous basin, both in terms of the 
nature of the ground (rather sorted sand grains) and of the vegetation. Except during 
flooding events, bed-load is the dominant mode of transport. The size of the grains on 
the river bed is around d — 330 /im. The experiments have been performed at the end 
of summer (low water period), in two straight and flat portions of the river. 

The experiments were conducted as follows. Using long parallel metallic bars, the 
surface of the sand bed was carefully leveled at time t = 0, and the formation of bedforms 
was directly observed. In order to measure the emergence of ripples, we used a (water 
proof) laser sheet inclined at a low angle to the horizontal. Taking pictures through a 
glass plate from the top, we have determined the height profile Z(x, t) along this line as a 
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Figure 14. Formation of ripples in a natural river, starting at t = from a flat sand bed. The 
experiment was performed in the Leyre river, at Sauniac bridge, on September 16 th 2008, for 
H — 52 cm, T = 0.21 and u, — 3 cm/s. The sediment is well sorted with a mean grain size 
d = 320 ±70 /im (a) Amplitude 2C, of the bed disturbances as function of time t. £ is determined 
from the auto-correlation of the bed profile Z(x,t). The solid line is the best fit by equation 
(5.2 1 and gives a = 3 10 -2 s _1 or equivalently a~ = 35 s. (b) Bed elevation profile measured 
by taking a picture of the bed enlightened by an inclined laser sheet. Low amplitude ripples can 
be detected from t — 25 s after flattening the sand bed. Non-linear effects make the amplitude 
saturate around t — 100 s and a clear avalanche slip face can be observed at t = 150 s. During 

90 mm does not evolve significantly (see Fig, 
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this linear stage, the wavelength A 

this initial phase, a pattern coarsening toward dunes is observed which saturates at A 
and £ = 4.5 mm after typically one hour: tracking the pattern during further ~ 14 hours, we 
did not observe any significant change of these characteristics. 



function of time. Figure |l4"p shows the evolution of one such a ripple from the initial stage 
where it is symmetric to the time at which an avalanche slip face develops. To determine 
the wavelength A and the amplitude 2C, of the ripple, we computed the auto-correlation of 
the profile C(8) = (Z(x)Z(x + 6)). Typically 25 s after the beginning of the experiment, 
C shows a secondary maximum whose position gives A and whose amplitude gives £. 
During the first 150 s, the wavelength does not evolve much whereas the amplitude 



grows and saturates (Fig. 15 d). As the final amplitude is not very large compared to the 
grain size d (and thus the initial noise level), it is difficult to get a convincing evidence of 
an exponential growth. Still, the curves we obtained are consistent with a linear regime 
over a factor of two in amplitude (Fig. |14[ ). From our non-linear analysis, we expect an 
amplitude equation of the form: 



whose solution is: 



dt 



c 



c 



Co 



y/l + exp(-2crt) 



(5.1) 



(5.2) 



One can see in figure 14 that the fit of this relation to the data is very good, so that 
observations are consistent with a formation of ripples by a linear instability saturated 
by non- linear hydrodynamical effects. Importantly, the rescaled growth rate a/(ku Sf ) is 
around 1CP 3 and the rescaled propagation speed w/(fcu*) around 10~ 2 . As shown in 
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Figure 15. (a-c) Formation of dunes in a natural river, starting at t = from a flat sand 
bed. The experiment was performed in the Leyre river, at Mios bridge, on September 17* 
2008, for H = 50 cm, T = 0.28 and u, = 4 cm/s. The grain size is d = 330 ± 70 /im. (a) 
Time evolution of the wavelength A. The pattern coarsening starts after 150 s and stops after 
~ 4000 s. (b) Same graph, but restricted to the linear regime (between t = and t — 150 s). 
(c) The photograph shows the dunes of wavelength 40 cm formed after 6000 s. (d) Formation of 
mega-dunes, starting from a flat sand bed. The experiment was performed in the Leyre river, 
at Mios bridge, for H — 44 cm, T = 0.30 and it, = 4 cm/s. The sand is polydisperse: it a 
mixture of sand grains of size ~ 330 ± 70 /im, which cover 60% of the surface, and of coarse 
grains larger than 600 /im, which represent 40% of the surface - but 9% of the grains and 60% 
of the mass. The photograph shows 3 m long mega-dunes with ~ 40 cm superimposed dunes, 
(e) Photograph of the Leyre river at Mios bridge showing the sharp transition between dunes 
(zone of medium sand) and mega-dunes (zone of medium and coarse sand mixed). 



Part 1, with such small dimensionless numbers, the motion of the bed can be ignored in 
the hydrodynamical treatment. 



28 A. Fourriere, P. Claudin and B. Andreotti 

5.2. Formation of dunes and mega-dunes 
We have observed the evolution of the patterns for typically two hours after flattening 



the sand bed (Fig. 15 i). A statistically stationary state is eventually reached, which 
corresponds to what was observed in the natural conditions, i.e. before the experi- 
ment (Fig. [15]:). As in fl ume experiments (|Venditti et al. 2005al IVenditti et al. 2005b[ 
Valance & Langlois 2005), we have observed a coarsening of the ripple pattern, i.e. a 
progressive growth of the wavelength by merging of bedforms (Raudkivi S z Witte 1 990. 



iRaud kivi 2006). Fi gure PHh shows that this growth is linear in time and stops when 
the wavelength A becomes on the order of the flow depth H . Both these processes 
and the time-scales over which they take place arc consistent with the observations of 
|Guy et al. 1966] for flume experiments at larger Froude numbers. Again, with a bed mo- 
tion at the scale of hours, the approximation of a flow over a steady relief is almost perfect . 
Consistently with the theory, we did not observe the emergence of wavelengths directly 
at the scale of the flow depttjf] We thus reach the conclusion that the formation of dunes 
should not be associated to a linear instability but to a non-linear pattern evolution. In 
the unbounded case (an infinite flow depth), it is probable that this pattern coarsen- 
ing would have no limit as it is driven by hydrodynamics, which is mostly self-similar. 
The observation that this coarsening stops at some final wavelength should therefore be 
associated to a stabilising mechanism, namely the presence of the free surface. 

The places selected for the experiments were of particular interest as different final 



wavelengths were observed across the river (Fig. 15:-e). The experiments just described 
above have been performed in of the side of the river (say, at a distance less than the 
third of the river width from the bank). In this 'external' region, the sand is well sorted 
and the dune wavelength is observed to saturate for a rescaled wave-number kH just 



below the resonant conditions (Fig. 16 1). By contrast, in the central part of the river, 



much larger bedforms are present (Fig. 151-e). They display superimposed dunes on 



their stoss side, and we call mega-dunes here-after. The river slope and the flow velocity 
were not significantly different in the dune and mega-dune regions. The major difference 
was the presence of coarse grains causing bed armouring in the central part of the river. 



As shown in figure 15 3, the transition between the two regions is rather sharp. We have 
flattened the bed over a zone of 12 m in length and 4 m in width to observe the formation 
of mega-dunes. The initial stage is the formation of ripples composed of small sand grains 
that merge, leading to the same dunes as described on the side of the river. However, 
in the course of this pattern coarsening, the inter-dune zone becomes richer in coarse 
grains so that the dunes eventually propagate on a bed which is more difficult to erode. 
They progressively amalgamate into mega-dunes of wavelength ten times larger than 
the dunes, covered with superimposed ripples and dunes. Even in the asymptotic state, 
superimposed bedforms are continuously generated. As they propagate faster than the 
mega-dune, they accumulate at its crest. During the transient of formation of mega-dunes 
(typically 5 hours in our experiments), the pattern is disordered and is not composed 
by a unique Fourier mode. As far as one can say without having explicitly performed 
a multi-scale analysis of the topography, structures of growing size were progressively 
formed, which become ordered as they reach the final mega-dune wavelength (between 
10-ff and 20H). 

In summary, a small difference in the experimental conditions (here, most probably, the 
presence of coarse grains or not) can significantly affect the pattern coarsening dynamics. 



t In the flume experiments of Carling et al. 2005] the large grain size (5 mm) and the mod- 
erate water depth (0.375 m) are such that the typical size of the first emerging bedforms (few 
decimeters) are comparable. 



29 




Figure 16. Representation in the T vs kH diagram of different data set. (a) Experiments per- 
formed in the Leyre river, starting from a fiat sand bed. (b) Flume experiments: wavelength of 
the stationary bedforms obtained at long time, (c) Bedforms in natural rivers. The mega-dunes 
are represented by white symbols whereas dunes are displayed with black symbols, (d) Water 
depth in rio Parana (H ~ 8 m, d ~ 300 |im, T = 0.16) measured by multibeam echo sounding 
(this photograph is from|Parsons et al. 2005). Mega-dunes of wavelength A ~ 125 m ~ 15-ff can 
be observed, with superimposed dunes of wavelength 6 m and probably superimposed ripples 
too (not visible). 



The wavelength of mega-dunes can be larger than those of dunes by one to two orders 



of magnitude (Fig. 16 i). With this observation, the question is not anymore to find 
a new destabilising mechanism to explain the dune formation. Rather, as larger and 
larger bedforms are produced by non-linear processes, one needs to identify a stabilising 
mechanism limiting this coarsening. As evidenced here, this is exactly the role played by 
the free surface. 

5.3. Conclusion 

Figure [16] presents several series of measurements of the final bedform wavelength col- 
lected in the literature for flume experiments (panel b) and natural rivers (panel c) . The 
most impressive data-set is certainly that reported by |Guy et al. 1966| As these au- 
thors separated their data into ripples and dunes on the basis of a definition that we do 
not agree with (see introduction), we have ignored here their annotations. Importantly, 
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|Guy et al. 1966| have not looked at the transient of formation of bedforms but have fo- 
cused on the long time regime. For example, they have often started a new experiment 
with a bed in the state reached at the end of the previous run. Moreover, the control 
parameters (slope and flow rate) were varied in the course of the experiment to main- 
tain constant secondary quantities such as the flow depth or the Froude number. This 
methodology is probably responsible for the lack of reproducibility and the huge disper- 
sion of data (one decade horizontally) on can see in figure 16 3. By contrast, the series 
of experiments performed by Baas 1999 are much more controlled and reproducible. In 
particular, the wavelength is measured as a function of time starting from a flat sand 
bed and fitted to obtain its asymptotic value. These flume experiments clearly evidence 
the difference between the initial wavelength at which ripples form, their evolution by 
pattern coarsening to form dunes and the non-linear selection of the final wavelength due 
to the free surface. Note that we have not taken into account the correction of raw mea- 
surements performed bv lBaas 19991 to take into account the temperature dependence of 
the viscosity, as it is a non-sense in the hydraulically rough regime. The last points in fig- 
ure 16 d have been obtained bv lRobert fc Uhlman 20011 As for |Guy et al. 1966 we have 
not taken into account the denomination of the bedforms (ripple or dunes) used in this 
article. The points obtained by Baas 19991 and IRobert fc Uhlman 2001 are very close to 
the resonant conditions; the whole data set of |Guy et al. 1966| is clearly in the subcritical 
regime - at low T and/or low wavenumber kH - and globally follows the resonance curve. 
We hav e gathered in figure [ 1 6^ the points measured in the Leyre river, in the Gradyd tidal 
chanel (Barthol dy fc al. 2005[ ), in the Mississipi river (IHarbor 1998jl . in the Missouri river 
(jAnnambhotla 1972j) , in the Rhine ( |Carling et al. 2000|IWilbers fc Ten Brinke 2003)) and 
in the Rio Parana (Par sons et al. 2005 jl . One can observe that the dunes propagating on 
the stoss slope of mega-dunes lie in the same region of the diagram as the simple dunes 
- roughly between kH = 0.1/JF 2 and kH = \/T 2 . The mega-dunes lie between kH = 0.3 
and kH = OA/ J 72 . It should be emphasised that the different points in this graph cor- 
respond to very different flow depth H (compare e.g. the megadunes in the Leyre river 
shown in figure [l4| i and those in the rio Parana shown in figure [l6ji) . It means that 
the saturation length, which determines the wavelength at which ripples form, is not a 
relevant length for the formation of dunes and mega-dunes. 

We then see why a distinction between ripples and dunes based on some absolute 
dimension considerations are misleading from the physical point of view: depending on 
the value of the water depth, bedforms 'of less than about 2 feet' ( |Guy et al. 1966| ), or 
'less than 0.6 m' ( |Ashley 1990| ), can feel or not the free surface. Nor a criterion based 
on the amplitude |<5|£ of distortion of the free surface can be satisfying to define dunes. 
Figure 17 shows the value of \5\ in the parameter space (T,kH), computed using the 
hydrodynamical model detailed in part 1. One observes that the free surface is distorted 
by less than 10% of the bedform amplitude for T < 0.35 or for kH > 3. With a 
definition based on the amplitude of distortion of the free surface such as that proposed 
by |Guy et al. 1966| there would be no dune at all in natural rivers! At low Froude 
numbers, dunes do not distort significantly the free surface. Approaching T — 1 from 
below, the distortion of the free surface, in antiphase with the topography, becomes more 
and more pronounced (Fig. 17 o). Above T ~ 0.6, this effect becomes sufficient to create 
a zone of stable wavelengths around the resonance (Fig. 13 1 and thus a gap to cross 
during the pattern coarsening. 

These results then suggest a new classification of river bedforms based on the dy- 
namical mechanisms responsible for their formation. The most obvious criterion is the 
sensitivity to the presence of the free surface. We thus define ripples as the bedforms 
whose wavelength A is sufficiently small compared to the flow depth H not to feel the 
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Figure 17. Amplitude \5\ of the distortion of the free surface related to that of the bedform. 
(a) Isocontours of |<J| in the plane T vs kH (color online), (b) Distortion |<S| computed along the 
resonance curve (black solid line in panel a) for H/zq = 10 2 (dashed line), H/zo = 10 3 (dotted 
line) and H/zq — 10 4 (solid line) as a function of the Froude number T. The amplitude in the 
resonant conditions is maximum at T ~ 1. 



finiteness of the flow depth. In the diagram T vs kH, they are located in the super- 
critical region, i.e. their wavenumber is larger than the resonant value. When the scale 
separation between the flow depth and the saturation length is sufficient, the structures 
that form by linear instability of a flat sand bed are ripples. Beyond the linear regime the 
pattern coarsening leads to growing wavelengths that reach the resonance curve. Once 
in the subcritical region on the left of this curve (see Fig. 16 1 hydrodynamics becomes 
affected and in some case dominated by the presence of the free surface and the bedforms 
can be called dunes. As evidenced by our field experiments, the pattern coarsening can 
end with very different bedforms depending on the conditions (e.g. the grain size distri- 
bution). The associated non-linear selection of pattern wavelength is an open problem 
(see lPoliti fc Misbah 2004llAndreotti fc al. 2006)) that will require specific investigations. 
Our observations suggest to define mega-dunes as bedforms sufficiently large to present 
superimposed dunes. With this definition, what is usually called 'river bars' in the litera- 
ture falls into our category of mega-dunes. Similarly, the term 'sand- wavelet' introduced 
by Coleman (see e.g. Coleman & Melville 1994, Coleman & Melville 1996) are ripples in 
their initial stage. It is furthermore worth noting that dunes can present superimposed 
ripples when the scale separation between the flow depth H and the saturation length 
-L sat is sufficient. The criterion to separate dunes from mega-dunes is then based on the 
location of the superimposed bedforms with respect to the resonance curve in the plane T 
vs kH. This distinction is important as in very deep water, the pattern coarsening would 
lead to ripples large enough to accommodate superimposed ripples on their stoss side. 
In this case, one should talk about mega-ripples. In summary, the appellation of given 
subaqueous bedforms should be chosen in function of their location in the diagram J- vs 
kH, with a suffix 'mega-' to point out the presence or not of superimposed structures. 



6. Summary of the results 

As the present paper aims both to present novel results and to discuss the state of the 
art of the physics of ripples and dunes, we wish to sum up in this section, in a qualitative 
manner, our methodology and key results. 
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We have shown that the linear instability of a flat sand bed can be abstracted in terms 
of four well-defined quantities: 

• the component of the basal shear stress in phase and in quadrature with the bottom 
(the parameters A and B, see part 1), 

• the relation between the sediment flux, the shear stress and its threshold value, which 
determines the time-scale of the instability but does not influence the initial wavelength, 

• the saturation length, i.e. the relaxation length between the sediment flux and the 
shear stress, which essentially determines the length-scale of the instability, 

• the relation between the threshold shear stress and the slope, which encodes the 
stabilising effect of gravity. 

The framework resulting from these ingredients is very general and can encompass a wide 
spectrum of specific models without any profound difference of nature. As discussed in 
the first part of these twin papers, the details of the hydrodynamical description can 
affect A and B to various degrees, depending for instance on the physics at work in the 
surface layer or on the way the different length-scales of the problem intermix with each 
other - see the summary section of part 1. From the choice of a particular transport 
model, one can compute the steady (or saturated) value of the sediment flux, as well as 
the characteristic length over which the actual flux relaxes to equilibrium. Finally, the 
expression of the threshold shear stress as a function of the grain size, the bed slope and 
other possible parameters (e.g. the cohesion between the grains) comes from a separate 
mechanical model. Doing so, one can then independently investigate the effect of each 
dynamical mechanism on the formation of ripples and dunes. 

We have devoted a large fraction of this part 2 to the study of two sediment transport 
regimes. They differ by the nature of the mechanisms that make the flux saturate: 

• an erosion limited regime, in which the negative feedback of the particle transport on 
the flow is negligible. The saturation of the flux is limited by the time needed to entrain 
one grain at the surface of the sand bed (jCharru 2006)) . The associated saturation length 
is determined by the distribution of the traps at the surface of the bed. It is proportional 
to the grain size and it should diverge when the shear velocity reaches the value for which 
the flow is sufficient to entrain all the grains at the surface. 

• a momentum limited regime, in which the flux saturates due to the extra-shear stress 
exerted by the moving grains on the fluid ( |Bagnold 19 56). If the shear velocity is not 
too large so that bed load is the dominant mode of transport, the saturation length is 
dominated by the grain inertia thus proportional to the density ratio between the grain 
and the fluid times the grain size. 

We have derived a detailed model of bed load describing both the erosion limited and the 

momentum limited regimes that is able to fit the data of Julicn 1998|and Fernandez Luque fc van Beek 1976 



We have also derived a quantitative model for the threshold Shields number, which takes 
into account the bed slope effect and is able to reproduce the observed decrease at the 
transition between the viscous and the turbulent regimes. Sheet flows, in which the 
grains flow over a significant number of grain diameters, and suspended load are outside 
the validity of the model because they require a different hydrodynamical treatement 
where the flow inside the transport layer must be modelled. 

These different ingredients at hand, we have performed the linear stability analysis of 
a flat sand bed. The destabilising mechanism is of hydrodynamical nature and is related 
to the phase advance of the basal shear stress with respect to the topography. Two 
stabilising mechanisms are identified: the sediment transport saturation length L sat and 
the slope effect, which depends on the ratio tt*/uth- As L sat is generically larger than 
zq, it dominates and essentially determines the scaling law followed by the most growing 
wavelength A max . This most unstable mode is associated to ripples, which thus form by 
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a linear instability. In the case of a smooth bottom, the roughness seen from the inner 
layer is governed by the viscous length v/u*, which may dominate the scaling of A max if 
vju* is larger than L sat QSumer fc Bakioglu 1984| . 

Due to the slope effect, the ratio A max /I/ sa t is a decreasing function of u*/uth- We 
have analysed different sets of measurements of initial ripple wavelengths A max available 
in the literature and deduced the corresponding saturation length, assuming that the hy- 
drodynamical model is correct. L sat is found roughly independent of and between 5 d 
and 15 d. These values are consistent with a saturation length limited by the grain iner- 
tia (L sa t — 2(Ps/pf) d) as previously stated by the authors (jClaudin fc Andreotti 2 006, 
lAndreott i & Claudin 2007). This prediction is especially good for large grains (d > 
1 mm), for which viscous effects are completely negligible. However, the data show sys- 
tematic dependencies that are not captured by the present model: the wavelength at 
which ripples form is systematically larger for smaller or rougher grains. 

We have performed field experiments in the Leyre river, whose results show that the 
evolution of the ripple wavelength A and amplitude £ at short time are consistent with a 
linear instability. We also observed that these ripples present a pattern coarsening: their 
wavelength grows and saturate just after crossing the resonance condition of surface 
waves. In the course of this pattern coarsening, the bedforms are in quasi-cquilibrium 
between erosion and deposition, a situation for which we have shown that the aspect 
ratio of ripples can be predicted quantitatively using the weakly non-linear description 
of hydrodynamics derived in part 1 . 

The influence of the river free surface on the bed is stabilising. As shown in part 1 and 
in agreement with standard hydraulics, the gravity waves excited at the free surface by 
the bedforms are in phase at small X/H (supercritical regime) and in antiphase at large 
X/H (subcritical regime). In between, the free surface is phase-advanced with respect to 
the bottom, so that it tends to induce a phase-delay of the shear stress on the ground. 
Bedforms of wavelength around the resonance conditions are thus stabilised. Moreover, 
at very large X/H, the inner layer invades the whole flow and the free surface again has 
a strong stabilising effect. As no destabilising mechanism is associated to the presence of 
the free surface, dunes do not form by a linear instability. This result is directly confirmed 
by our field experiments showing that they form by non-linear pattern coarsening, as 
suggested by Raudkivi & Wittc 1990 and Raudkivi 2006 Finally, our experiments show 
that the non-linear selection of the final wavelength is very sensitive to small changes in 
the experimental conditions, and in particular to the presence of coarse grains. 

Finally, our results suggest to classify subaqueous bedforms according to the dynam- 
ical mechanisms that control their formation. We thus propose the following bedform 
definitions and characteristics: 

• Ripples are bedforms whose wavelength A is sufficiently small compared to the flow 
depth H not to feel the presence of the free surface. In the diagram T vs kH, they 
are located on the right of the resonance curve (supercritical regime). They form by 
linear instability and their initial wavelength essentially scales on the saturation length 
L sat . They exhibit pattern coarsening and remain ripples until they cross the condition 
of resonance of the surface waves. 

• Dunes are bedforms whose wavelength A is sufficiently large compared to the flow 
depth H to be stabilised by the presence of the free surface. The ripple pattern coarsening 
generically leads to the formation of dunes. In the diagram T vs kH, they are located 
along the resonance curve on the subcritical side. If the flow depth H is much larger 
than the wavelength at which ripples form, dunes may present superimposed ripples. 

• Like dunes, mega-dunes are under the influence of the free surface and on the left 
side of the resonance curve but they present superimposed dunes. They typically result 
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from the coarsening of a dune pattern pushed to very large wavelength by heterogeneities 
(in particular a polydispersed sediment). 

Further studies arc needed to investigate this non-linear wavelength selection in details. 
In the purpose of describing the interactions between bedforms, future models will have to 
incorporate hydrodynamical non-linearities and in particular flow separation, for which 
the systematic expansion with respect to the amplitude performed in part 1 may serve 
as a well-controlled starting point. 



The field work has been carried out in agreement with the Pare Naturel Regional des 
Landes de Gascogne. This work has benefited from the financial support of the french 
minister of research. 



Appendix A. Static threshold 

In this appendix, we derive an explicit expression of the transport threshold, valid 
both in the viscous and turbulent regimes. We introduce the effective velocity U of the 
flow around the sand grains and define for convenience its normalised counterpart: 

U= - = (Al) 

v(ps/pf - l )9 d 

For the sake of simplicity, we will take for U the fluid velocity at the altitude z = d/2. A 
good approximation of the relation between the shear velocity and the typical velocity 
around the grain U is given by 

ul = ^U+ 2 , , U 2 . (A2) 
d ln 2 (l + l/2r) V ; 

where v is the kinematic viscosity of the fluid, k the von Karman constant and r the 
aerodynamic roughness rescaled by the grain diameter d. The rescaled bed roughness r 
is on the order of r = 1/10. 

At the threshold, the horizontal force balance on a grain of the bed reads 

^( Ps - Pf )gd 3 = *-C dPf Uld\ (A3) 

where /i is the avalanche slope for sand grains, i.e. tan 32°. Cd is the drag coefficient, 
which is a function of the grain Reynolds number TZ. With a good accuracy, the drag 
law for natural grains can be put under the form: 

C d = (CU 2 + sll- 1 ' 2 ) 2 with K=™. (A 4) 

At this stage, we introduce the viscous size d v , defined as the diameter of grains whose 
free fall Reynolds number Uf a \\d/v is unity: 

d v = {p a /p I -\)-V* v 2 l*g-'/\ (A5) 

It corresponds to a grain size at which viscous and gravity effects are of the same order 
of magnitude. From the three previous relations, we get the equation for UtW- 

C^ th + S (^) 3/4 <^(f) 1/2 , (A6) 
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Figure 18. Geometry of the trajectory calculation. 



which solves into 
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The corresponding expression of the static threshold Shields number is finally: 



e th = 2 



d v 



3/2 



Uv 



d J ' 1^(1 + 1/2?-)" 

In the viscous regime, the above relations simplify into: 
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,3s 2 J \d„ 
In the turbulent regime, we get: 
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and 6 t h = 



4/iK 2 



3C 0o ln' ! (H-l/2r) 



(A 7) 



(A 8) 



(A 9) 



(A 10) 



As can be seen in figure[2] the most interesting range of grain diameters, between 100 /im 
and 1 mm, is precisely the zone of transition between the two regimes. The crossover 
grain diameter d co scales as: 

4o= ozbF^' (A11) 

which is around 200 /im. 



Appendix B. Motion of a single grain 

To describe the elementary jump by one grain diameter, we write the equation of 
motion of one grain dragged by a flow of effective velocity U and submitted to gravity. 
We assume that it looses completely its energy during collisions, due to the thin layer 
of fluid between the grains. Contrarily to Bagnold 1956 and Charm ct al. 2008, we do 
not introduce any 'effective friction': instead, we take into account the geometry of the 



trapping grains (Fig. 18 1, which is the physical effect responsible for this friction, as for 



the transport threshold. 

Starting from the position of the grain dug, we get its velocity V — d8u e+7r / 2 and 
its tangential acceleration d9ug + ^^ 2 - We assume that the reaction of the substrate is 
normal to it. Using, as in the appendix [X] d as a unit size and [(p s /pf — ^gd] 1 ^ 2 as a 
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unit velocity (sec Eq. Al), the dynamical equation reads: 



Cli 2 [9 2 -2cos9U9 + U 2 



1/4 



3/4 



(cos9U - 9j + sinf 



(Bl) 



2 Go ° 



(9 2 - 2cos9U9+U 2 



i/4 



1/2 
l.h 



+ 3- 1/ X 1/2 



(cos(9W-fl) +miW4. 



where [i has been taken equal to tan(7r/6) = l/VS. As in the appendix |a| the rescaled 
velocity U is related to the Shields number by: 



e = 2 'i 



3/2 



ln 2 (l + l/2r) 



(B2) 



This equation can be integrated numerically between 9 = — 7r/6 and = 7r/6, starting 
from = 0. Note that there is no adjustable parameter in this description. 

In order to get the asymptotic behavior analytically, we perform an expansion with 
respect to 9 around — ir/6 and with respect to U around Uth- We get a linear equation 
of the form: 



1 



o5/4 

+ K^CU 2 \9- 



tii 



2 / 7r\ 2n(U-U th ) 



with 



VsV'' 6 

8tt 



&/ 2 aU th 



(B3) 



3 3 / 2 (2 + 3^C^ 2 U th ) ' 

The solution of this linear equation is: 

7T Tr(U th -U) [ r_exp(r + t) -r + cxp(r_£) 

(7 = 1 1 H 

6 3aZYth 

where r_ and r + are the solutions of the equation: 



r + — r_ 



0. 



The time T needed to pass over d thus scales around the threshold as: 

uUth 



T In 



1 



1 



U-U th 



1 Pfd 
(ps - Pf)g 



(B4) 



(B5) 



(B6) 



Making use of equation (A 2 1, one sees that this erosion time diverges at the threshold 
as: 

Tex -In [6 -e th ]. (B7) 

On the other hand, T should resume to d/U far from the threshold. Integrating numer- 
ically equation (B 1 1 , we have found convenient approximations that follow the correct 
asymptotic behaviors just above the threshold and far from it. In the turbulent regime, 
it reads: 



T cx In 



/Q-ye 



th 



th 



' Pfd 
(ps - Pf)g 



(Bl 



e 


+ e th " 




.© 


-e th . 






37 



(B9) 



where the pref actor is around 1.5. In the viscous regime, it reads 

3/2 

T oc In 
with a prefactor around 17. 



Appendix C. Saturation length in a continuum description of 
momentum limited bed load 

We start from the momentum balance equation, as written bv lParker 19751 

% + J k) dx< *) = PfU * ~ Tsand( - 9 ^ ( C ^ 

By definition, the flux q saturates to g sa t when the shear stress is equal at the top of the 
transport layer and below it, on the bottom: 

Tsand(<7sat(w*)) = Pfu\. (C 2) 

Linearising the dynamical equation, we get: 



d, Q'sat ^ \ ^sand 
tq+ -r—o x q = 



feat-<z). (C3) 

« = 9aat 



h J dq 

Identifying this equation with the definition of the saturation length L sat and the satu- 
ration time T sat , 

T sat d t q + L sat d x q = q sat - <7, (C 4) 

we get: 

T ysat Ps Qsat 

ho (it pf2h u* du» 
One readily see that L sa t, as predicted by this expression, vanishes at the transport 
threshold as g sat and increases very fast with u*. None of these behaviors is correct, 
whatever the shape of the bed load equation q sa t(u*) ~ Einstein-Brown, Bagnold or 
those derived in the text. There are two simple reasons for this failure. First, the 
continuum approach ignores the process by which the momentum is transferred from 
the fluid to the grains (namely the drag of the grains). Second, the flux q evolves both 
because the number of moving particles evolves and because their velocity changes. A 
valid description should thus both present an equation for the erosion and a description 
of the trajectories, as done in the text. 
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